﻿using System;
using System.Text;
using SoftwareConsulting.BI11.Studio.InternalServices;

namespace SoftwareConsulting.BI11.NumericalMethods.MatrixAlgebra
{
    /// <summary>
    /// Класс декомпозиции нижней-верхней перестановки
    /// </summary>
    public class LowerUpperPermutationDecomposition
    {        
        private double[,] _rows;     
        private int[] _permutation = null;        
        private int _parity = 1;//четность перестановки

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="components"></param>
        public LowerUpperPermutationDecomposition(double[,] components)
        {
            if (components.GetLength(0) != components.GetLength(1))
                throw new BIMatrixNotSquareException(components.GetLength(0), components.GetLength(1), null);
            _rows = components;
            Initialize();
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="matrix"></param>
        public LowerUpperPermutationDecomposition(Matrix matrix)
        {
            if (!matrix.IsSquare)
                throw new BIMatrixNotSquareException(matrix.RowsCount, matrix.ColumnsCount, null);
            Initialize(matrix.Components);
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="matrix"></param>
        public LowerUpperPermutationDecomposition(SymmetricMatrix matrix)
        {
            Initialize(matrix.Components);
        }

        /// <summary>
        /// Возвращает обратную подстановку
        /// </summary>
        /// <param name="xTilde"></param>
        /// <returns></returns>
        private double[] GetBackwardSubstitution(double[] xTilde)
        {
            int n = _rows.GetLength(0);
            double[] result = new double[n];
            for (int i = n - 1; i >= 0; i--)
            {
                result[i] = xTilde[i];
                for (int j = i + 1; j < n; j++)
                    result[i] -= _rows[i, j] * result[j];
                result[i] /= _rows[i, i];
            }
            return result;
        }

        /// <summary>
        /// Выполняет декомпозицию
        /// </summary>
        private void Decompose()
        {
            int n = _rows.GetLength(0);
            _permutation = new int[n];
            for (int i = 0; i < n; i++)
                _permutation[i] = i;
            _parity = 1;
            try
            {
                for (int i = 0; i < n; i++)
                {
                    SwapRows(i, GetLargestPivotIndex(i));
                    Pivot(i);
                }
            }
            catch (ArithmeticException ex) 
            {
                Trace.StartBlock("Exception killed");
                Trace.Message("Arithmetic exception cathed while lower-upper permutation decomposition, and killed. Exception text:");
                Trace.Message(ex.ToString());
                Trace.EndBlock();
                _parity = 0; 
            } ;
        }

        /// <summary>
        /// Возвращает флаг проведенной декомпозиции
        /// </summary>
        /// <returns></returns>
        private bool IsDecomposed
        {
            get
            {
                if (_parity == 1 && _permutation == null)
                    Decompose();
                return _parity != 0;
            }
        }

        /// <summary>
        /// Возвращает детерминант матрицы
        /// </summary>
        /// <returns></returns>
        public double GetDeterminant()
        {
            if (!IsDecomposed)
                return double.NaN;
            double determinant = _parity;
            for (int i = 0; i < _rows.GetLength(0); i++)
                determinant *= _rows[i, i];
            return determinant;
        }

        /// <summary>
        /// Возвращает прямую подстановку
        /// </summary>
        /// <param name="c"></param>
        /// <returns></returns>
        private double[] GetForwardSubstitution(double[] c)
        {
            int n = _rows.GetLength(0);
            double[] result = new double[n];
            for (int i = 0; i < n; i++)
            {
                result[i] = c[_permutation[i]];
                for (int j = 0; j <= i - 1; j++)
                    result[i] -= _rows[i, j] * result[j];
            }
            return result;
        }

        /// <summary>
        /// Выполняет инициализацию вычислений
        /// </summary>
        private void Initialize()
        {
            _permutation = null;
            _parity = 1;
        }

        /// <summary>
        /// Выполняет инициализацию вычислений
        /// </summary>
        private void Initialize(double[,] components)
        {
            _rows = (double[,])components.Clone();
            Initialize();
        }

        /// <summary>
        /// Возвращает инверсированные компоненты матрицы
        /// </summary>
        /// <returns></returns>
        public double[,] GetInversedMatrixComponents()
        {
            if (!IsDecomposed)
                return null;
            int n = _rows.GetLength(0);
            double[,] inverseRows = new double[n, n];
            double[] column = new double[n];
            for (int i = 0; i < n; i++)
            {
                Array.Clear(column, 0, n);
                column[i] = 1;
                column = GetSolution(column);
                for (int j = 0; j < n; j++)
                    inverseRows[i, j] = column[j];
            }
            return inverseRows;
        }

        /// <summary>
        /// Возвращает индекс наибольшей оси вращения
        /// </summary>
        /// <param name="columnIndex"></param>
        /// <returns></returns>
        private int GetLargestPivotIndex(int columnIndex)
        {
            double maximum = Math.Abs(_rows[columnIndex, columnIndex]);
            double abs;
            int index = columnIndex;
            for (int i = columnIndex + 1; i < _rows.GetLength(0); i++)
            {
                abs = Math.Abs(_rows[i, columnIndex]);
                if (abs > maximum)
                {
                    maximum = abs;
                    index = i;
                }
            }
            return index;
        }

        /// <summary>
        /// Вращает относительно элемента с указанным индексом
        /// </summary>
        /// <param name="columnIndex"></param>
        private void Pivot(int columnIndex)
        {
            double inversePivot = 1 / _rows[columnIndex, columnIndex];
            int columnIndex2 = columnIndex + 1;
            int n = _rows.GetLength(0);
            for (int i = columnIndex2; i < n; i++)
            {
                _rows[i, columnIndex] *= inversePivot;
                for (int j = columnIndex2; j < n; j++)
                    _rows[i, j] -= _rows[i, columnIndex] * _rows[columnIndex, j];
            }
        }

        /// <summary>
        /// Возвращает решение
        /// </summary>
        /// <param name="c"></param>
        /// <returns></returns>
        public double[] GetSolution(double[] c)
        {
            return IsDecomposed
                            ? GetBackwardSubstitution(GetForwardSubstitution(c))
                            : null;
        }

        /// <summary>
        /// Возвращает решение
        /// </summary>
        /// <param name="cVector"></param>
        /// <returns></returns>
        public Vector GetSolution(Vector cVector)
        {
            double[] components = GetSolution(cVector.Components);
            return components == null ? null : new Vector(components);
        }

        /// <summary>
        /// Меняет i-ю и k-ю строки местами
        /// </summary>
        /// <param name="i"></param>
        /// <param name="k"></param>
        private void SwapRows(int i, int k)
        {
            if (i == k)
                return;            
            for (int j = 0; j < _rows.GetLength(1); j++)
            {
                double temp = _rows[i, j];
                _rows[i, j] = _rows[k, j];
                _rows[k, j] = temp;
            }
            int nTemp = _permutation[i];
            _permutation[i] = _permutation[k];
            _permutation[k] = nTemp;
            _parity = -_parity;            
        }

        /// <summary>
        /// Делает компоненты матрицы симметричными
        /// </summary>
        /// <param name="components"></param>
        public static void SymmetrizeComponents(double[,] components)
        {
            for (int i = 0; i < components.GetLength(0); i++)
            {
                for (int j = i + 1; j < components.GetLength(1); j++)
                {
                    components[i,j] += components[j,i];
                    components[i,j] *= 0.5;
                    components[j,i] = components[i,j];
                }
            }
        }

        /// <summary>
        /// Возвращает как строку
        /// </summary>
        /// <returns></returns>
        public override string ToString()
        {
            StringBuilder sb = new StringBuilder();
            char[] separator = { '[', ' ' };
            int n = _rows.GetLength(0);
            for (int i = 0; i < n; i++)
            {
                separator[0] = '{';
                for (int j = 0; j < n; j++)
                {
                    sb.Append(separator);
                    sb.Append(_rows[i, j]);
                    if (j == 0) separator[0] = ' ';
                }
                sb.Append('}');
                sb.Append('\n');
            }
            if (_permutation != null)
            {
                sb.Append(_parity == 1 ? '+' : '-');
                sb.Append("( " + _permutation[0]);
                for (int i = 1; i < n; i++)
                    sb.Append(", " + _permutation[i]);
                sb.Append(')');
                sb.Append('\n');
            }
            return sb.ToString();
        }
    }
}
